Multinomial Flight Modeling

Flight Trials Winter 2020 were conducted from 2/17/2020 - 3/10/2020. Soapberry bugs were flight tested twice (T1 vs. T2) for multiple hours in the flight mill and observed from 8 AM to (5-8 PM) each day. I used multicategorical, nominal modeling (multinom) to analyze how sex and host plant as well as changes in mass and egg-laying changed a soapberry bug’s flight response between trials.

Delta Flight Response (T1 vs. T2)

Introduction

Clean the Data

Flight Case ~ Mass Diff

Flight Case ~ Mass Per

Flight Case ~ Mass Diff or Mass Perc + Sex

Flight Case ~ Mass Diff or Mass Perc + Sex + Host Plant

Flight Case ~ Mass Diff or Mass Perc + Sex + Wing2Body

Females Only

Introduction

Graphs I generated below illustrate how changes in mass (more specifically % mass) or egg-laying led to changes in flight response. To generate these graphs, I ran multicategorical regressions to model the probability of different flight delta cases due to sex, host plant, and changes in mass or egg-laying. Changes in speed and distance as a result of changes in mass or egg-laying were not analyzed here but rather in a separate script.

To perform these analyses, a new variable will be created called “flight_case” which calculates the flight response differences between T1 and T2. Here is a formula and key describing each flight case event and encoding:

Formula : flight case = response in \(T_2\) - response in \(T_1\)

Delta Flight Response Key
Event Encoding
flew in both trials 2
flew in T2 only 1
flew in neither trials 0
flew in T1 only -1

Since the outcomes (or response variables) were no longer binomial, I used multicategorical logit models (refer to Chapter 6 of An Introduction to Categorical Data Analysis, Second Edition by Alan Agresti). Below I’ve written out notes on performing these analyses.

What we are interested in is a multicategorical, nominal response variable. When the response variable is nominal, there is no natural order among the response variable categories (unordered categories). When the response variable is multicategorical, its multicategory models assume that the counts in the categories of \(Y\) have a multinomial distribution.

Let J = number of categories for \(Y\).

\(\pi_1,...\pi_J\) = the response probabilities where \(\sum_j \pi_j = 1\).

With \(n\) independent observations, the probability distribution for the number of outcomes of the J types is the multinomial. It specifies the probability for each possible way the \(n\) observations can fall in the J categories. Multicategory logit models simultaneously use all pairs of categories by specifying the odds of outcome in one category instead of another.

Logit models for nominal response variables pair each category with a baseline category. The choice of the baseline category is arbitrary. When the last category (J) is the baseline, the baseline-category logits are

\(\log (\frac{\pi_j}{\pi_J}), j = 1, ..., J - 1\).


Given that the response falls in category j or category J, this is the log odds that the response is j. For J = 3, for instance, the model uses \(\log(\pi_1/\pi_3)\) and \(\log(\pi_2/\pi_3)\). The baseline-category logit model with a predictor x is

\(\log \frac{\pi_j}{\pi_J} = \alpha_J + \beta_j x, j = 1, ..., J - 1\)


The model has J - 1 equations, with separate parameters for each. The effects vary according to the category paired with the baseline. When J = 2, this model simplifies to a single equation for \(\log(\pi_1/\pi_2) = logit(\pi_1)\), resulting in ordinary logistic regression for binary responses.

So how do these pair of categories determine equations for all other pairs of categories? Here is an arbitrary pair of categories a and b that follows the general equation above,

\(\log (\frac{\pi_a}{\pi_b}) = \log ( \frac{\pi_a/pi_J}{\pi_b/pi_J}) = \log (\frac{\pi_a}{\pi_J}) - \log (\frac{\pi_b}{\pi_J})\)

\(= (\alpha_a + \beta_a x) - (\alpha_b + \beta_b x)\)

\(= (\alpha_a - \alpha_b) + (\beta_a - \beta_b) x\)


So, the equation for categories a and b has the form \(\alpha + \beta x\) with intercept parameter \(\alpha = (\alpha_a - \alpha_b)\) and with slope parameter \(\beta = (\beta_a - \beta_b)\).

Let’s apply it to our data now:

Cleaning the Data

source_path = "~/Desktop/git_repositories/SBB-dispersal/avbernat_working_on/Rsrc/"

script_names = c("center_flight_data.R", # Re-centers data 
                 "clean_flight_data.R", # Loads and cleans data
                 "unique_flight_data.R",
                 "get_warnings.R", 
                 "compare_models.R",
                 "regression_output.R", # Cleans regression outputs; prints in color or B&W
                 "AICprobabilities.R")

for (script in script_names) { 
  path = paste0(source_path, script)
  source(path) 
}
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
output_col = FALSE # Change to TRUE if working in Base R or RStudio; FALSE if generating an HTML

data <- read_flight_data("data/all_flight_data-Winter2020.csv")
data_all <- data[[1]]
data_tested <- data[[2]]
d <- create_delta_data(data_tested)
## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
colnames(d)[c(1:2,5,66:68,71:73)]
## [1] "ID"          "sex"         "host_plant"  "egg_diff"    "mass_diff"  
## [6] "flew_diff"   "mass_per"    "flight_case" "egg_case"

Age Effect

Age also plays a role in this dataset, but age was unknown because the soapberry bugs were field collected.

Modeling

Below I used the multinom function from the nnet package to estimate a multinomial logistic regression model. There are other functions in other R packages capable of multinomial regression. I chose the multinom function because it does not require the data to be reshaped (as the mlogit package does) and to mirror the example code found in Hilbe’s Logistic Regression Models.

First, I chose the baseline by specifying the level using relevel function. Then, I ran our model using multinom. The multinom package does not include p-value calculations for the regression coefficients, so I calculated P using Wald tests (i.e. z-tests).

Finally, I considered mass difference as my predictor variable where a (+) sign means the bug gained mass between trials and a (-) sign means the bug lost mass between trials.

Delta Mass Key
Event Sign
gained mass from T1 to T2 pos
no mass change between trails 0
lost mass from T1 to T2 neg

Baseline

df <- d %>%
  filter(!is.na(mass_diff), !is.na(flight_case)) 
  
df <- df[with(df, order(mass_diff)),]
n_trials = nrow(df)

df$flight_case <- relevel(as.factor(df$flight_case), ref = "0")

Null model

null <- multinom(flight_case ~ 1, data = df) 
## # weights:  8 (3 variable)
## initial  value 385.389832 
## iter  10 value 319.269929
## final  value 319.269680 
## converged

Multinom model: flight_case ~ mass_diff

model <- multinom(flight_case ~ mass_diff, data = df) 
## # weights:  12 (6 variable)
## initial  value 385.389832 
## iter  10 value 317.510452
## iter  20 value 310.451788
## iter  30 value 308.978871
## iter  40 value 308.544941
## iter  50 value 308.493264
## iter  60 value 308.481503
## iter  70 value 308.459043
## iter  80 value 308.445007
## iter  90 value 308.440565
## final  value 308.439954 
## converged
model_table = calculate_P(model)
## 
##  AIC:  628.8799 
##    (Intercept)  Estimate DF Std. Err.     z   wald P > |z|
## -1      -0.770    69.211  6    16.478 4.200 17.641   0.000
## 1       -2.129     6.699  6    28.958 0.231  0.054   0.817
## 2        0.420    25.328  6    12.662 2.000  4.002   0.045
anova(null, model, test="Chisq") # adding mass_diff improves fit 
## Likelihood ratio tests of Multinomial Models
## 
## Response: flight_case
##       Model Resid. df Resid. Dev   Test    Df LR stat.      Pr(Chi)
## 1         1       831   638.5394                                   
## 2 mass_diff       828   616.8799 1 vs 2     3 21.65945 7.678753e-05

Which multinomial model equations are significant:

get_significant_models()
## [1] -1
##   (Intercept)  Estimate DF Std. Err.         z      wald      P > |z|
## 0   0.7621587 -67.22835  6  16.32315 -4.118588 16.962770 3.812004e-05
## 1  -1.3799779 -54.68207  6  30.71245 -1.780453  3.170012 7.500189e-02
## 2   1.1808369 -42.94895  6  13.95307 -3.078100  9.474699 2.083251e-03
## [1] 0
##    (Intercept)  Estimate DF Std. Err.         z        wald      P > |z|
## -1  -0.7701329 69.210777  6  16.47825 4.2001301 17.64109245 2.667617e-05
## 1   -2.1286461  6.699274  6  28.95809 0.2313438  0.05351994 8.170477e-01
## 2    0.4195114 25.328344  6  12.66161 2.0004041  4.00161667 4.545664e-02
## [1] 1
##    (Intercept)  Estimate DF Std. Err.          z      wald    P > |z|
## 0     2.134069 -12.69279  6  29.55112 -0.4295196 0.1844871 0.66754511
## -1    1.366080  56.16348  6  30.76233  1.8257221 3.3332614 0.06789213
## 2     2.555220  12.15032  6  29.18084  0.4163802 0.1733724 0.67713184
## [1] 2
##    (Intercept)  Estimate DF Std. Err.          z       wald     P > |z|
## 1    -2.549995 -17.17484  6  28.83335 -0.5956589  0.3548096 0.551403105
## 0    -0.420632 -24.93984  6  12.65319 -1.9710315  3.8849650 0.048720277
## -1   -1.194381  44.61978  6  14.07959  3.1691113 10.0432666 0.001529058

Estimated Parameters, Standard Errors, and Wald Test Statistics for All Main Effects Model
Odds Effect Parameter Estimate SE exp(\(\frac{\beta}{50}\)) Wald p
\(\frac{P(Yi=flew\,twice)}{P(Yi=flew\, in\, T1\, only)}\) Intercept \(\alpha\) 1.18
\(\delta\) Mass (g) \(\beta_1\) -42.95 13.95 0.42 9.47 0
\(\frac{P(Yi=flew\,twice)}{P(Yi=did\, not\, fly)}\) Intercept \(\alpha\) 0.42
\(\delta\) Mass (g) \(\beta_1\) 25.33 12.66 1.66 0.17 0.05

The multinomial (ML) prediction equations are:

prediction_equations(model_table) 
## [1] "log(pi_-1 / pi_1) = 1.36 + 62.51 Mass Change     Flew in T1, rather than T2"   
## [2] "log(pi_2 / pi_-1) = 1.19 + -43.88 Mass Change     Flew in both, rather than T1"
## [3] "log(pi_2 / pi_1) = 2.55 + 18.63 Mass Change     Flew in both, rather than T2"
N Model P
278 Flew in T1 only rather than T2 only = 1.36 + 62.51 Mass Change not sig
278 Flew in both rather than T1 only = 1.19 - 43.88 Mass Change sig
278 Flew in both rather than T2 only = 2.55 + 18.63 Mass Change not sig

Estimated Odds

For bugs of mass_diff x + 1/50 (0.02) g, the estimated odds that the bug flew twice rather than only in T1 equals

exp(-43.88/50)
## [1] 0.4157796

times the estimated odds at mass_diff x (g). As gain mass, then less likely to fly again.

exp(coef(model)/50) # this compares to not flying, the baseline | 1/0 not significant
##    (Intercept) mass_diff
## -1   0.9847154  3.991693
## 1    0.9583206  1.143376
## 2    1.0084255  1.659584

Predicted probabilities

You can also use predicted probabilities to help you understand the model. The multicategory logit model has an alternative expression in terms of the response probabilities. This is

\(\pi_j = \frac{e^{\alpha_j + \beta_j x}}{\sum_he^{\alpha_h + \beta_h x}}, j=1,...,J\)


The denominator is the same for each probability, and the numerators for various j sum to the denominator where \(\sum_j \pi_j = 1\). The parameters (\(\alpha\) and\(\beta\)) equal zero for whichever category is the baseline in the logit expressions. So these would be the equations for the model above:

\(\hat{\pi_{-1}} = \frac{e^{-1.77 + 69.21 x}}{1 + e^{-1.77 + 69.21 x} + e^{-2.13 - 6.70 x} + e^{0.42 + 25.33 x}}, j=1,...,J\)

\(\hat{\pi_1} = \frac{e^{-2.13 - 6.70 x}}{1 + e^{-1.77 + 69.21 x} + e^{-2.13 - 6.70 x} + e^{0.42 + 25.33 x}}, j=1,...,J\)

\(\hat{\pi_2} = \frac{e^{0.42 + 25.33 x}}{1 + e^{-1.77 + 69.21 x} + e^{-2.13 - 6.70 x} + e^{0.42 + 25.33 x}}, j=1,...,J\)

\(\hat{\pi_0} = \frac{e^{0 + 0 x} = 1}{1 + e^{-1.77 + 69.21 x} + e^{-2.13 - 6.70 x} + e^{0.42 + 25.33 x}}, j=1,...,J\)

x = mass_diff, which can give you the estimated probabilities.


You can calculate these predicted probabilities for each of our outcome levels using the fitted function. You can start by generating the predicted probabilities for the observations in our dataset and viewing the first few rows. Then, you can plot them.

head(pp <- fitted(model))
##           0         -1          1         2
## 1 0.6268782 0.01288591 0.05518211 0.3050538
## 2 0.5919733 0.01843237 0.05424679 0.3353476
## 3 0.5859452 0.01955212 0.05405531 0.3404474
## 4 0.5612480 0.02470149 0.05318314 0.3608674
## 5 0.5485567 0.02772717 0.05268168 0.3710345
## 6 0.5356441 0.03109397 0.05213548 0.3811264

Multinom model: flight case ~ mass_per

Re-order

df <- df[with(df, order(mass_per)),]
model <- multinom(flight_case ~ mass_per, data = df) 
## # weights:  12 (6 variable)
## initial  value 385.389832 
## iter  10 value 307.053645
## final  value 306.984778 
## converged
model_table = calculate_P(model)
## 
##  AIC:  625.9696 
##    (Intercept)  Estimate DF Std. Err.     z   wald P > |z|
## -1      -0.909     0.044  6     0.010 4.492 20.179   0.000
## 1       -2.132     0.001  6     0.020 0.043  0.002   0.965
## 2        0.344     0.021  6     0.008 2.494  6.221   0.013

Model comparing flying only in T2 to the baseline (not flying at all) is not significant (P = 0.97). That is probably because so few bugs (only male) flew in T2 only.

Significant ML equations

run_multinom_model = function(d) {
  m <- multinom(flight_case ~ mass_per, trace=FALSE, data = d) 
  return(m)
}
get_significant_models()
## [1] -1
##   (Intercept)    Estimate DF   Std. Err.         z      wald      P > |z|
## 0   0.9091574 -0.04408961  6 0.009815063 -4.492035 20.178382 7.054571e-06
## 1  -1.2221618 -0.04330188  6 0.020160403 -2.147868  4.613337 3.172424e-02
## 2   1.2530401 -0.02354676  6 0.007814186 -3.013336  9.080191 2.583931e-03
## [1] 0
##    (Intercept)     Estimate DF   Std. Err.          z         wald      P > |z|
## -1  -0.9093175 0.0440898971  6 0.009814924 4.49212803 20.179214198 7.051501e-06
## 1   -2.1317451 0.0008628504  6 0.019875088 0.04341366  0.001884746 9.653718e-01
## 2    0.3436623 0.0205503655  6 0.008238958 2.49429177  6.221491411 1.262088e-02
## [1] 1
##    (Intercept)      Estimate DF  Std. Err.           z        wald    P > |z|
## 0     2.127991 -0.0008619011  6 0.01984339 -0.04343517 0.001886614 0.96535464
## -1    1.217852  0.0432625817  6 0.02012413  2.14978612 4.621580358 0.03157214
## 2     2.471342  0.0196994818  6 0.01941135  1.01484321 1.029906739 0.31018057
## [1] 2
##    (Intercept)    Estimate DF   Std. Err.         z     wald     P > |z|
## 1   -2.4737239 -0.01975061  6 0.019437023 -1.016134 1.032527 0.309565810
## 0   -0.3436951 -0.02053882  6 0.008239362 -2.492768 6.213892 0.012675166
## -1  -1.2532127  0.02356475  6 0.007815213  3.015241 9.091678 0.002567752

Estimated Parameters, Standard Errors, and Wald Test Statistics for All Main Effects Model
Odds Effect Parameter Estimate SE exp(\(\beta*20\)) Wald p
\(\frac{P(Yi=flew\,twice)}{P(Yi=flew\, in\, T1\, only)}\) Intercept \(\alpha\) 1.25
\(\delta\) Mass % \(\beta_1\) -0.02 0.01 0.62 9.08 0
\(\frac{P(Yi=flew\,twice)}{P(Yi=did\, not\, fly)}\) Intercept \(\alpha\) 0.34
\(\delta\) Mass % \(\beta_1\) 0.02 0.01 1.51 6.22 0.01
\(\frac{P(Yi=flew\, in\, T1\, only)}{P(Yi=flew\, in\, T2\, only)}\) Intercept \(\alpha\) 1.22
\(\delta\) Mass % \(\beta_1\) 0.04 0.02 2.38 4.62 0.03

The multinomial (ML) prediction equations are:

gsub("Mass Change",  "Mass Percent Change", prediction_equations(model_table))
## [1] "log(pi_-1 / pi_1) = 1.22 + 0.04 Mass Percent Change     Flew in T1, rather than T2"   
## [2] "log(pi_2 / pi_-1) = 1.25 + -0.02 Mass Percent Change     Flew in both, rather than T1"
## [3] "log(pi_2 / pi_1) = 2.48 + 0.02 Mass Percent Change     Flew in both, rather than T2"
N Model P
278 Flew in T1 only rather than T2 only = 1.22 + 0.04 Mass % sig
278 Flew in both rather than T1 only = 1.25 - 0.02 Mass % sig
278 Flew in both rather than T2 only = 2.48 + 0.02 Mass % not sig

Estimated Odds

For bugs of mass_perc x + 20%*(0.04) the estimated odds that bug flew in T1 rather than T2 equals

exp(0.04*20)
## [1] 2.225541

times the estimated odds at mass_perc x (%). So for every 20% increase in a bug’s original mass, the bug is 2.23 times more likely to fly only in T1. That makes sense because the odds of flying in T2 only are very low to begin with.

For bugs of mass_perc x + 40%*(-0.02), the estimated odds that the bug flew twice instead of only in T1 equals

exp(40*-0.02) # 2.5 times less likely to fly twice if gain mass
## [1] 0.449329
exp(-40*-0.02) # 2.5 times more likely to fly twice if loose mass
## [1] 2.225541

times the estimated odds at mass_perc x (%). Similar message - for every 40% increase in a bug’s original mass, the bug was 2.23 times more likely to fly only in T1 rather than fly twice. Although, when a bug looses 40% of original mass then becomes more likely to fly twice then only once.

exp(coef(model)*20) # this compares to no fly, the baseline # gain large % mass
##     (Intercept) mass_per
## -1 1.264064e-08 2.415238
## 1  3.047172e-19 1.017407
## 2  9.660792e+02 1.508336
exp(coef(model)*5) # gain small % mass
##     (Intercept) mass_per
## -1 1.060333e-02 1.246637
## 1  2.349493e-05 1.004324
## 2  5.575107e+00 1.108216
exp(coef(model)*-5) # loose small % mass
##     (Intercept)  mass_per
## -1 9.431000e+01 0.8021582
## 1  4.256237e+04 0.9956950
## 2  1.793688e-01 0.9023509
exp(coef(model)*-20) # loose large % mass
##     (Intercept)  mass_per
## -1 7.910994e+07 0.4140378
## 1  3.281732e+18 0.9828910
## 2  1.035112e-03 0.6629821
# 1/0 not significant

1 vs. 0 is not significant as already stated. However, -1 vs. 0 shows that as gain more mass then become increasingly more likely to only fly once. If loose mass then most likely to not fly at all, and more likely to fly twice then fly only once.

head(pp <- fitted(model))
##           0         -1          1         2
## 1 0.5386483 0.04284419 0.06190317 0.3566043
## 2 0.5356362 0.04375052 0.06158899 0.3590243
## 3 0.5290568 0.04577594 0.06090110 0.3642662
## 4 0.5263205 0.04663701 0.06061436 0.3664282
## 5 0.4982068 0.05614792 0.05764729 0.3879980
## 6 0.4949167 0.05734314 0.05729764 0.3904425

For more on multinomial plotting for when have more than 1 predictor variables see the following: https://stats.idre.ucla.edu/r/dae/multinomial-logistic-regression/

Now let’s compare some models

First: flight_case, mass_diff, and sex

df <- df[with(df, order(mass_diff)),]
data <- data.frame(R = df$flight_case, 
         A = df$mass_diff,
         B = df$sex_c)
source("src/compare_models.R")
model_comparisonsAIC("src/generic multinomial models- multinom 1RF + 2 FF.R")
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
##        [,1]      [,2]     
## AICs   589.7091  593.7062 
## models 3         4        
## probs  0.8804523 0.1193253
## 
## m3   multinom(formula = R ~ A + B, data = data, trace = FALSE)
## m4   multinom(formula = R ~ A * B, data = data, trace = FALSE)
anova(m3, m4, test="Chisq") # Adding A*B does not improve fit
## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##   Model Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1 A + B       825   571.7091                                
## 2 A * B       822   569.7062 1 vs 2     3 2.002834 0.5718186
delta_mass_model <- multinom(flight_case ~ mass_diff + sex_c, data = df)
## # weights:  16 (9 variable)
## initial  value 385.389832 
## iter  10 value 297.481645
## iter  20 value 287.067261
## iter  30 value 286.119187
## iter  40 value 285.989229
## iter  50 value 285.976242
## iter  60 value 285.937258
## iter  70 value 285.917216
## iter  80 value 285.877940
## iter  80 value 285.877937
## iter  90 value 285.874179
## iter 100 value 285.854534
## final  value 285.854534 
## stopped after 100 iterations
model_table = calculate_P2(delta_mass_model, "mass_diff", "sex_c")
## 
##  AIC:  589.7091 
##    (Intercept) mass_diff   sex_c DF    SE1   SE2     z1      z2 ML wald
## -1      -0.965    70.907  -0.788  9 17.662 0.219  4.015  -3.590  16.117
## 1      -17.461   -27.383 -16.279  9 67.108 0.184 -0.408 -88.332   0.167
## 2        0.205    20.417  -0.885  9 12.907 0.157  1.582  -5.626   2.502
##    sex wald P1 > |z| P2 > |z|
## -1   12.886    0.000        0
## 1  7802.495    0.683        0
## 2    31.647    0.114        0

Significant models

MASS_ML = get_significant_models(11) # mass_diff

MASS_ML
## [[1]]
##   (Intercept) mass_diff        sex_c DF      SE1       SE2         z1
## 0   0.9417796 -66.75437   0.78783396  9 17.26359 0.2181963 -3.8667728
## 1 -16.7482516 -62.19962 -15.68367725  9 64.88462 0.1921608 -0.9586187
## 2   1.1562140 -47.99979  -0.09980787  9 16.89219 0.2153509 -2.8415381
##            z2    ML wald     sex wald     P1 > |z|     P2 > |z|
## 0   3.6106663 14.9519321   13.0369114 0.0001102851 0.0003054114
## 1 -81.6174659  0.9189499 6661.4107356 0.3377508448 0.0000000000
## 2  -0.4634662  8.0743389    0.2148009 0.0044896483 0.6430302553
## 
## [[2]]
##    (Intercept) mass_diff       sex_c DF      SE1       SE2         z1
## -1  -0.9654638  70.90708  -0.7875466  9 17.66235 0.2193886  4.0145884
## 1  -17.4607727 -27.38313 -16.2794390  9 67.10782 0.1842989 -0.4080468
## 2    0.2045214  20.41677  -0.8851817  9 12.90706 0.1573493  1.5818294
##            z2    ML wald   sex wald     P1 > |z|     P2 > |z|
## -1  -3.589734 16.1169202   12.88619 5.954962e-05 3.310159e-04
## 1  -88.331734  0.1665022 7802.49529 6.832393e-01 0.000000e+00
## 2   -5.625585  2.5021842   31.64721 1.136885e-01 1.848803e-08
## 
## [[3]]
##    (Intercept)  mass_diff    sex_c DF      SE1       SE2         z1       z2
## 0     14.87658 -10.791307 13.64842  9 64.16531 0.2084871 -0.1681798 65.46410
## -1    13.91935  60.297896 12.84533  9 64.69628 0.2241167  0.9320149 57.31539
## 2     15.08554   9.494493 12.75626  9 63.59561 0.1997793  0.1492948 63.85179
##       ML wald sex wald  P1 > |z| P2 > |z|
## 0  0.02828443 4285.549 0.8664419        0
## -1 0.86865178 3285.053 0.3513288        0
## 2  0.02228893 4077.051 0.8813210        0
## 
## [[4]]
##    (Intercept) mass_diff      sex_c DF      SE1        SE2         z1
## 1   -7.1936471 -34.99467 -4.8981646  9 65.80680 40.1975619 -0.5317789
## 0   -0.2077583 -17.72600  0.8831577  9 12.71691  0.1568779 -1.3938920
## -1  -1.1597837  50.20531  0.1029897  9 17.04736  0.2155423  2.9450486
##            z2   ML wald    sex wald    P1 > |z|     P2 > |z|
## 1  -0.1218523 0.2827888  0.01484798 0.594879142 9.030160e-01
## 0   5.6295860 1.9429348 31.69223846 0.163350233 1.806427e-08
## -1  0.4778168 8.6733114  0.22830892 0.003229039 6.327806e-01
SEX_ML = get_significant_models(12) # sex

SEX_ML
## [[1]]
##   (Intercept) mass_diff        sex_c DF      SE1       SE2         z1
## 0   0.9417796 -66.75437   0.78783396  9 17.26359 0.2181963 -3.8667728
## 1 -16.7482516 -62.19962 -15.68367725  9 64.88462 0.1921608 -0.9586187
## 2   1.1562140 -47.99979  -0.09980787  9 16.89219 0.2153509 -2.8415381
##            z2    ML wald     sex wald     P1 > |z|     P2 > |z|
## 0   3.6106663 14.9519321   13.0369114 0.0001102851 0.0003054114
## 1 -81.6174659  0.9189499 6661.4107356 0.3377508448 0.0000000000
## 2  -0.4634662  8.0743389    0.2148009 0.0044896483 0.6430302553
## 
## [[2]]
##    (Intercept) mass_diff       sex_c DF      SE1       SE2         z1
## -1  -0.9654638  70.90708  -0.7875466  9 17.66235 0.2193886  4.0145884
## 1  -17.4607727 -27.38313 -16.2794390  9 67.10782 0.1842989 -0.4080468
## 2    0.2045214  20.41677  -0.8851817  9 12.90706 0.1573493  1.5818294
##            z2    ML wald   sex wald     P1 > |z|     P2 > |z|
## -1  -3.589734 16.1169202   12.88619 5.954962e-05 3.310159e-04
## 1  -88.331734  0.1665022 7802.49529 6.832393e-01 0.000000e+00
## 2   -5.625585  2.5021842   31.64721 1.136885e-01 1.848803e-08
## 
## [[3]]
##    (Intercept)  mass_diff    sex_c DF      SE1       SE2         z1       z2
## 0     14.87658 -10.791307 13.64842  9 64.16531 0.2084871 -0.1681798 65.46410
## -1    13.91935  60.297896 12.84533  9 64.69628 0.2241167  0.9320149 57.31539
## 2     15.08554   9.494493 12.75626  9 63.59561 0.1997793  0.1492948 63.85179
##       ML wald sex wald  P1 > |z| P2 > |z|
## 0  0.02828443 4285.549 0.8664419        0
## -1 0.86865178 3285.053 0.3513288        0
## 2  0.02228893 4077.051 0.8813210        0
## 
## [[4]]
##    (Intercept) mass_diff      sex_c DF      SE1        SE2         z1
## 1   -7.1936471 -34.99467 -4.8981646  9 65.80680 40.1975619 -0.5317789
## 0   -0.2077583 -17.72600  0.8831577  9 12.71691  0.1568779 -1.3938920
## -1  -1.1597837  50.20531  0.1029897  9 17.04736  0.2155423  2.9450486
##            z2   ML wald    sex wald    P1 > |z|     P2 > |z|
## 1  -0.1218523 0.2827888  0.01484798 0.594879142 9.030160e-01
## 0   5.6295860 1.9429348 31.69223846 0.163350233 1.806427e-08
## -1  0.4778168 8.6733114  0.22830892 0.003229039 6.327806e-01
prediction_equations2(model_table, " Mass Change", " Sex ")
## [1] "Where F = 1"
## [1] "log(pi_-1 / pi_1) = 16.5 + 98.29 Mass Change + 15.49 Sex      Flew in T1, rather than T2"  
## [2] "log(pi_2 / pi_-1) = 1.17 + -50.49 Mass Change + -0.1 Sex      Flew in both, rather than T1"
## [3] "log(pi_2 / pi_1) = 17.67 + 47.8 Mass Change + 15.39 Sex      Flew in both, rather than T2"
N Model P
278 Flew in T1 only rather than T2 only = 16.5 + 98.29 Mass Change + 15.49 Sex not mass | sex**
278 Flew in both rather than T1 only = 1.25 - 50.49 Mass Change - 0.1 Sex mass** | not sex
278 Flew in both rather than T2 only = 17.67 + 47.8 Mass Change + 15.39 Sex not mass | sex**

Second: flight_case, mass_per, and sex

df <- df[with(df, order(mass_per)),]
data <- data.frame(R = df$flight_case, 
         A = df$mass_per,
         B = df$sex_c)
source("src/compare_models.R")
model_comparisonsAIC("src/generic multinomial models- multinom 1RF + 2 FF.R")
##        [,1]      [,2]      
## AICs   587.5607  593.4209  
## models 3         4         
## probs  0.9492355 0.05068255
## 
## m3   multinom(formula = R ~ A + B, data = data, trace = FALSE)
## m4   multinom(formula = R ~ A * B, data = data, trace = FALSE)
anova(m3, m4, test="Chisq") # Adding A*B does not improve fit
## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##   Model Resid. df Resid. Dev   Test    Df  LR stat.   Pr(Chi)
## 1 A + B       825   569.5607                                 
## 2 A * B       822   569.4209 1 vs 2     3 0.1398496 0.9866598
delta_mass_model <- multinom(flight_case ~ mass_per + sex_c, data = df)
## # weights:  16 (9 variable)
## initial  value 385.389832 
## iter  10 value 286.869825
## iter  20 value 284.809036
## iter  30 value 284.797822
## final  value 284.780360 
## converged
model_table = calculate_P2(delta_mass_model, "mass_per", "sex_c")
## 
##  AIC:  587.5607 
##    (Intercept) mass_per  sex_c DF   SE1   SE2     z1      z2 ML wald sex wald
## -1      -1.015    0.043 -0.692  9 0.010 0.203  4.390  -3.408  19.272   11.617
## 1       -6.820   -0.009 -5.626  9 0.026 0.183 -0.348 -30.721   0.121  943.764
## 2        0.124    0.019 -0.902  9 0.008 0.159  2.334  -5.684   5.447   32.310
##    P1 > |z| P2 > |z|
## -1    0.000    0.001
## 1     0.728    0.000
## 2     0.020    0.000
prediction_equations2(model_table, " Mass Percent Change", " Sex ")
## [1] "Where F = 1"
## [1] "log(pi_-1 / pi_1) = 5.81 + 0.05 Mass Percent Change + 4.93 Sex      Flew in T1, rather than T2"    
## [2] "log(pi_2 / pi_-1) = 1.14 + -0.02 Mass Percent Change + -0.21 Sex      Flew in both, rather than T1"
## [3] "log(pi_2 / pi_1) = 6.94 + 0.03 Mass Percent Change + 4.72 Sex      Flew in both, rather than T2"

Significant models

run_multinom_model = function(d) {
  m <- multinom(flight_case ~ mass_per + sex_c, trace=FALSE, data = d) 
  model_table = calculate_P2(m, "mass_per", "sex_c", print_table=FALSE)
  return(model_table)
}
MASS_ML = get_significant_models(11) # mass per

MASS_ML
## [[1]]
##   (Intercept)    mass_per      sex_c DF         SE1       SE2        z1
## 0    1.015581 -0.04252208  0.6918554  9 0.009686166 0.2029804 -4.389981
## 1   -6.371413 -0.05173676 -5.5003717  9 0.026394139 0.1894281 -1.960161
## 2    1.139438 -0.02314856 -0.2096559  9 0.008374107 0.1965823 -2.764303
##           z2   ML wald   sex wald     P1 > |z|     P2 > |z|
## 0   3.408483 19.271930  11.617760 1.133607e-05 0.0006532504
## 1 -29.036729  3.842231 843.131617 4.997699e-02 0.0000000000
## 2  -1.066504  7.641369   1.137432 5.704460e-03 0.2861957009
## 
## [[2]]
##    (Intercept)     mass_per      sex_c DF         SE1       SE2         z1
## -1  -1.0154279  0.042521764 -0.6917921  9 0.009686178 0.2029708  4.3899423
## 1   -6.8204622 -0.009162539 -5.6256871  9 0.026356158 0.1831234 -0.3476432
## 2    0.1236475  0.019385661 -0.9015721  9 0.008305979 0.1586103  2.3339406
##            z2    ML wald  sex wald     P1 > |z|     P2 > |z|
## -1  -3.408333 19.2715935  11.61673 1.133807e-05 6.536117e-04
## 1  -30.720742  0.1208558 943.76396 7.281082e-01 0.000000e+00
## 2   -5.684196  5.4472786  32.31009 1.959883e-02 1.314291e-08
## 
## [[3]]
##    (Intercept)    mass_per    sex_c DF        SE1       SE2        z1       z2
## 0     6.872767 0.009217925 5.678090  9 0.02636245 0.1992609 0.3496612 28.49576
## -1    5.857132 0.051742281 4.986267  9 0.02639275 0.2103403 1.9604732 23.70572
## 2     6.996531 0.028594653 4.776587  9 0.02573154 0.1901576 1.1112685 25.11910
##      ML wald sex wald  P1 > |z| P2 > |z|
## 0  0.1222629 812.0081 0.7265930        0
## -1 3.8434553 561.9610 0.0499405        0
## 2  1.2349177 630.9692 0.2664528        0
## 
## [[4]]
##    (Intercept)    mass_per      sex_c DF         SE1       SE2        z1
## 1   -6.9585308 -0.02867121 -4.7375707  9 0.025752878 0.1697712 -1.113321
## 0   -0.1237005 -0.01938216  0.9023844  9 0.008306705 0.1586318 -2.333316
## -1  -1.1391262  0.02314430  0.2098559  9 0.008374251 0.1965973  2.763745
##            z2  ML wald   sex wald    P1 > |z|     P2 > |z|
## 1  -27.905622 1.239483 778.723750 0.265570750 0.000000e+00
## 0    5.688546 5.444362  32.359553 0.019631589 1.281258e-08
## -1   1.067441 7.638288   1.139429 0.005714213 2.857729e-01
SEX_ML = get_significant_models(12) # sex

SEX_ML
## [[1]]
##   (Intercept)    mass_per      sex_c DF         SE1       SE2        z1
## 0    1.015581 -0.04252208  0.6918554  9 0.009686166 0.2029804 -4.389981
## 1   -6.371413 -0.05173676 -5.5003717  9 0.026394139 0.1894281 -1.960161
## 2    1.139438 -0.02314856 -0.2096559  9 0.008374107 0.1965823 -2.764303
##           z2   ML wald   sex wald     P1 > |z|     P2 > |z|
## 0   3.408483 19.271930  11.617760 1.133607e-05 0.0006532504
## 1 -29.036729  3.842231 843.131617 4.997699e-02 0.0000000000
## 2  -1.066504  7.641369   1.137432 5.704460e-03 0.2861957009
## 
## [[2]]
##    (Intercept)     mass_per      sex_c DF         SE1       SE2         z1
## -1  -1.0154279  0.042521764 -0.6917921  9 0.009686178 0.2029708  4.3899423
## 1   -6.8204622 -0.009162539 -5.6256871  9 0.026356158 0.1831234 -0.3476432
## 2    0.1236475  0.019385661 -0.9015721  9 0.008305979 0.1586103  2.3339406
##            z2    ML wald  sex wald     P1 > |z|     P2 > |z|
## -1  -3.408333 19.2715935  11.61673 1.133807e-05 6.536117e-04
## 1  -30.720742  0.1208558 943.76396 7.281082e-01 0.000000e+00
## 2   -5.684196  5.4472786  32.31009 1.959883e-02 1.314291e-08
## 
## [[3]]
##    (Intercept)    mass_per    sex_c DF        SE1       SE2        z1       z2
## 0     6.872767 0.009217925 5.678090  9 0.02636245 0.1992609 0.3496612 28.49576
## -1    5.857132 0.051742281 4.986267  9 0.02639275 0.2103403 1.9604732 23.70572
## 2     6.996531 0.028594653 4.776587  9 0.02573154 0.1901576 1.1112685 25.11910
##      ML wald sex wald  P1 > |z| P2 > |z|
## 0  0.1222629 812.0081 0.7265930        0
## -1 3.8434553 561.9610 0.0499405        0
## 2  1.2349177 630.9692 0.2664528        0
## 
## [[4]]
##    (Intercept)    mass_per      sex_c DF         SE1       SE2        z1
## 1   -6.9585308 -0.02867121 -4.7375707  9 0.025752878 0.1697712 -1.113321
## 0   -0.1237005 -0.01938216  0.9023844  9 0.008306705 0.1586318 -2.333316
## -1  -1.1391262  0.02314430  0.2098559  9 0.008374251 0.1965973  2.763745
##            z2  ML wald   sex wald    P1 > |z|     P2 > |z|
## 1  -27.905622 1.239483 778.723750 0.265570750 0.000000e+00
## 0    5.688546 5.444362  32.359553 0.019631589 1.281258e-08
## -1   1.067441 7.638288   1.139429 0.005714213 2.857729e-01
Estimated Parameters, Standard Errors, and Wald Test Statistics for All Main Effects Model
Odds Effect Parameter Estimate SE exp(\(\beta\)) Wald p
\(\frac{P(Yi=flew\,twice)}{P(Yi=flew\, in\, T1\, only)}\) Intercept \(\alpha\) 1.14
\(\delta\) Mass % \(\beta_1\) -0.02 0.01 0.98 7.64 0.01
Sex \(\beta_2=F\) -0.21 0.2 0.81 1.14 0.29
\(\frac{P(Yi=flew\,twice)}{P(Yi=flew\, in\, T2\, only)}\) Intercept \(\alpha\) 7
\(\delta\) Mass % \(\beta_1\) 0.03 0.03 1.03 1.23 0.27
Sex \(\beta_2=F\) 4.78 0.19 118.7 630.97 0
\(\frac{P(Yi=flew\,twice)}{P(Yi=did\, not\, fly)}\) Intercept \(\alpha\) 0.12
\(\delta\) Mass % \(\beta_1\) 0.02 0.01 1.02 5.45 0.02
Sex \(\beta_2=F\) -0.9 0.16 0.41 32.31 0
\(\frac{P(Yi=flew\, in\, T1\, only)}{P(Yi=flew\, in\, T2\, only)}\) Intercept \(\alpha\) 5.86
\(\delta\) Mass % \(\beta_1\) 0.05 0.03 46.69 4.62 0.05
Sex \(\beta_2=F\) 4.99 0.21 146.39 561.96 0
\(\frac{P(Yi=flew\, in\, T1\, only)}{P(Yi=did\, not\, fly)}\) Intercept \(\alpha\) -1.02
\(\delta\) Mass % \(\beta_1\) 0.04 0.01 1.04 19.27 0
Sex \(\beta_2=F\) -0.69 0.2 0.5 11.62 0
\(\frac{P(Yi=flew\, in\, T1\, only)}{P(Yi=did\, not\, fly)}\) Intercept \(\alpha\) -6.82
\(\delta\) Mass % \(\beta_1\) -0.01 0.03 0.99 0.12 0.73
Sex \(\beta_2=F\) -5.63 0.18 0 943.76 0
N Model P
278 Flew in T1 only rather than T2 only = 5.81 + 0.05 Mass % + 4.93 Sex mass** | sex**
278 Flew in both rather than T1 only = 1.14 - 0.02 Mass % - 0.21 Sex mass** | not sex
278 Flew in both rather than T2 only = 6.94 + 0.03 Mass % + 4.72 Sex not mass | sex**

Estimated Odds

For males:

For bugs of mass_perc x + 20%*(0.05), the estimated odds that bug flew in T1 rather than T2 equals

exp(0.05*20)
## [1] 2.718282

times the estimated odds at mass_perc x (%). So for every 20% increase in a bug’s original mass, the bug was 2.72 times more likely to fly only in T1. That makes sense because the odds of flying in T2 only are very low to begin with.

For bugs of mass_perc x + 20*(-0.02%), the estimated odds that the bug flew twice instead of only in T1 equals

exp(20*-0.02) # 20 times less likely to fly twice if gain mass
## [1] 0.67032
exp(-20*-0.02) # 20 times more likely to fly twice if loose mass
## [1] 1.491825

times the estimated odds at mass_perc x (%). Similar message - for every 20% increase in a bug’s original mass, the bug was 1.5 times more likely to fly only in T1. So, when a bug starts loosing mass it becomes more likely to fly twice then only once.

For females:

For bugs of mass_perc x + 20%*(0.05), the estimated odds that bug flew in T1 rather than T2 equals

exp(0.05*20 + 4.93)
## [1] 376.1545

times the estimated odds at mass_perc x (%). So for every 20% increase in a bug’s original mass, the bug was 376 times more likely to fly only in T1. Crazy! It doesn’t take much for a female to only fly once.

For bugs of mass_perc x + 20%*(0.05), the estimated odds that bug flew in both trials rather than T1 equals

exp(-0.02*20 - 0.21)
## [1] 0.5433509
exp(-0.02*-40 - 0.21)
## [1] 1.803988

times the estimated odds at mass_perc x (%). So half as likely to fly in both trials if gain mass. Would need to loose a lot of mass (40) to be almost twice as likely (1.8) to fly in both trials.

exp(coef(delta_mass_model)*20) # this compares to no fly, the baseline # gain large % mass
##     (Intercept)  mass_per        sex_c
## -1 1.513930e-09 2.3406655 9.798741e-07
## 1  5.730838e-60 0.8325593 1.367423e-49
## 2  1.185738e+01 1.4736071 1.475858e-08
exp(coef(delta_mass_model)*5) # gain small % mass
##     (Intercept)  mass_per        sex_c
## -1 6.237728e-03 1.2369007 3.146245e-02
## 1  1.547229e-15 0.9552209 6.081010e-13
## 2  1.855655e+00 1.1017814 1.102202e-02
exp(coef(delta_mass_model)*-5) # loose small % mass
##     (Intercept)  mass_per        sex_c
## -1 1.603148e+02 0.8084723 3.178392e+01
## 1  6.463168e+14 1.0468783 1.644464e+12
## 2  5.388934e-01 0.9076211 9.072748e+01
exp(coef(delta_mass_model)*-20) # loose large % mass
##     (Intercept)  mass_per        sex_c
## -1 6.605326e+08 0.4272289 1.020539e+06
## 1  1.744946e+59 1.2011156 7.313027e+48
## 2  8.433568e-02 0.6786069 6.775721e+07
# 1/0 mass not sig but sex**
  • If F and gain a lot of mass then most likely to fly only in T1 rather than none.
  • If F and loose mass then most likely to not fly in either trial.

Predicted Probabilities

head(pp <- fitted(delta_mass_model))
##           0         -1            1         2
## 1 0.7917303 0.03003973 4.362037e-06 0.1782256
## 2 0.7894639 0.03073036 4.325625e-06 0.1798015
## 3 0.7844677 0.03228066 4.247094e-06 0.1832474
## 4 0.7823713 0.03294258 4.214827e-06 0.1846819
## 5 0.7601763 0.04036342 3.895616e-06 0.1994564
## 6 0.7574981 0.04130977 3.859619e-06 0.2011882
df$index = 1:nrow(df)
females = df %>%
  filter(sex=="F")
males = df %>%
  filter(sex=="M")
Frows = females$index
Mrows = males$index

No females flew in T2 only.

Adding Host Plant

df <- df[with(df, order(mass_diff)),]
data <- data.frame(R = df$flight_case, 
         A = df$mass_diff,
         B = df$sex_c,
         C = df$host_c)
source("src/compare_models.R")
model_comparisonsAIC("src/generic multinomial models- multinom 1RF + 3 FF.R")
##        [,1]      [,2]      [,3]       [,4]       [,5]     
## AICs   589.7091  591.9755  593.0435   593.7062   593.9087 
## models 4         12        13         8          7        
## probs  0.5142205 0.1655767 0.09707198 0.06969088 0.0629815
## 
## m4   multinom(formula = R ~ A + B, data = data, trace = FALSE)
## m12  multinom(formula = R ~ A * C + B, data = data, trace = FALSE)
## m13  multinom(formula = R ~ B * C + A, data = data, trace = FALSE)
## m8   multinom(formula = R ~ A * B, data = data, trace = FALSE)
## m7   multinom(formula = R ~ A + B + C, data = data, trace = FALSE)
anova(m4, m7, test="Chisq") # Adding C (host plant) does not improve fit
## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##       Model Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1     A + B       825   571.7091                                
## 2 A + B + C       822   569.9087 1 vs 2     3 1.800378 0.6148527

Host plant is not significant.

Wing2Body

df <- df[with(df, order(mass_per)),]
df$wing2body_c = df$wing2body - mean(df$wing2body)
data <- data.frame(R = df$flight_case, 
         A = df$mass_per,
         B = df$sex_c,
         C = df$wing2body_c)
model_script = paste0(source_path,"generic multinomial models- multinom 1RF + 3 FF.R")
model_comparisonsAIC(model_script)
##        [,1]      [,2]      [,3]      
## AICs   582.2678  585.1197  587.133   
## models 7         12        13        
## probs  0.6671688 0.1603139 0.05858546
## 
## m7   multinom(formula = R ~ A + B + C, data = data, trace = FALSE)
## m12  multinom(formula = R ~ A * C + B, data = data, trace = FALSE)
## m13  multinom(formula = R ~ B * C + A, data = data, trace = FALSE)
anova(m7, m12, test="Chisq") # adding A*C does not improve fit
## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##       Model Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1 A + B + C       822   558.2678                                
## 2 A * C + B       819   555.1197 1 vs 2     3 3.148182 0.3693379
anova(m7, m13, test="Chisq") # Adding B*C does not improve fit
## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##       Model Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1 A + B + C       822   558.2678                                
## 2 B * C + A       819   557.1330 1 vs 2     3 1.134887 0.7686596
calculate_P3 = function(m, print_table=TRUE) {
  
  # compute p-values using Wald tests (i.e. z-tests)
  s <- summary(m) 
  z <- s$coefficients/s$standard.errors 
  wald <- z^2
  p <- (1 - pnorm(abs(z), 0, 1)) * 2
    
  # summary table (similar to lm, glm, or lmer function)
  model_table <- cbind(s$coefficients, c(s$edf, s$edf, s$edf), 
                       s$standard.errors[,1:4], z[,1:4], wald[,1:4], p[,1:4])
  colnames(model_table) <- c("(Intercept)", "mass %", "sex","wing2body","DF", "SEi",
                              "SE1", "SE2", "SE3","zi", "z1", "z2", "z3",
                              "waldi", "wald1","wald2", "wald3", "Pi>|z|", "P1>|z|", "P2>|z|", "P3>|z|")
  if (print_table){
    cat("\n", "AIC: ", s$AIC, "\n") 
    print(round(model_table,3))    
  }
  return(model_table)
}
#model <- multinom(flight_case ~ mass_per + sex_c + wing2body, data = df)
#model_table = calculate_P3(model)
model <- multinom(flight_case ~ mass_per + sex_c + wing2body_c, data = df)
## # weights:  20 (12 variable)
## initial  value 385.389832 
## iter  10 value 286.740091
## iter  20 value 280.436850
## iter  30 value 279.437125
## iter  40 value 279.174660
## iter  50 value 279.134087
## final  value 279.133921 
## converged
model_table = calculate_P3(model)
## 
##  AIC:  582.2678 
##    (Intercept) mass %    sex wing2body DF   SEi   SE1   SE2    SE3      zi
## -1      -0.935  0.041 -0.571    23.739 12 0.243 0.010 0.212 12.059  -3.854
## 1       -8.177 -0.005 -6.954    -6.595 12 0.187 0.025 0.187 18.786 -43.767
## 2        0.201  0.018 -0.760    28.094 12 0.172 0.008 0.166  9.718   1.173
##        z1      z2     z3    waldi  wald1    wald2 wald3 Pi>|z| P1>|z| P2>|z|
## -1  4.254  -2.698  1.969   14.850 18.096    7.278 3.875  0.000  0.000  0.007
## 1  -0.215 -37.102 -0.351 1915.510  0.046 1376.550 0.123  0.000  0.830  0.000
## 2   2.141  -4.590  2.891    1.375  4.585   21.071 8.357  0.241  0.032  0.000
##    P3>|z|
## -1  0.049
## 1   0.726
## 2   0.004

Signal of flying in only T2 is weak, which is probably why the equation in the model was not significant. It does not differ much from not flying at all.

Significant equations

run_multinom_model = function(d) {
  m <- multinom(flight_case ~ mass_per + sex_c + wing2body_c, trace=FALSE, data = d)
  model_table = calculate_P3(m, print_table=FALSE)
  return(model_table)
}

# ML's below are all the same
MASS_PER_ML = get_significant_models(19) # mass%

SEX_ML = get_significant_models(20) # sex

WING2BODY_ML = get_significant_models(21) # wing2body

getting_odds = function(Odds, effect, pars, odds, r) {
  # Odds:     Equation of the odds written out as vector of strings.
  # effect:   Main effects written out as a vector of strings.
  # pars:     Parameters written out as a vector of strings
  # odds:     Model table ML equations for the given outcome and baseline.
  # r:        Row in which the outcome is indexed. 
  estimate = round(unlist(odds[r,1:4]),2)
  SE = c(round(unlist(odds[r,6:9]),2))
  for (i in 2:4) {
    if (i == 2) {
      estimate[i] = estimate[i]*20
    }
    if (i == 4) {
      estimate[i] = round(estimate[i]/100,2)
    }
  }
  exp.b = c(round(exp(estimate[1:4]),2))
  wald = c(round(unlist(odds[r,14:17]),2))
  pvals = c(round(unlist(odds[r,18:21]),3))
  
  for (i in 1:length(pvals)){
    if (pvals[i] == "0") {
      pvals[i] = "<0.001"
    }
  }
  key = cbind(Odds, effect, pars, estimate, SE, exp.b, wald, pvals) 
  
  return(key)
}
# 2 vs. -1
#Odds = c("$\\frac{P(Yi=flew\\,twice)}{P(Yi=flew\\, in\\, T1\\, only)}$", " ", " ", " ")
Odds = c("$\\underline{P(Yi=flew\\,twice)}$", "$P(Yi=flew\\, in\\, T1\\, only)$", " ", " ")
effect = c("Intercept", "$\\delta$ Mass %", "Sex", "Wing-to-Body")
pars = c("$\\alpha$", "$\\beta_1$", "$\\beta_2\\, (F=1 | M=-1)$", "$\\beta_3\\,(\\mu=0)$ ")
key1 = getting_odds(Odds, effect, pars, MASS_PER_ML[[1]], 3) 

# 2 vs. 1
#Odds = c("$\\frac{P(Yi=flew\\,twice)}{P(Yi=flew\\, in\\, T2\\, only)}$", " ", " ", " ")
Odds = c("$\\underline{P(Yi=flew\\,twice)}$", "$P(Yi=flew\\, in\\, T2\\, only)$", " ", " ")
key2 = getting_odds(Odds, effect, pars, MASS_PER_ML[[3]], 3) 

# 2 vs. 0
#Odds = c("$\\frac{P(Yi=flew\\,twice)}{P(Yi=did\\, not\\, fly)}$", " ", " ", " ")
Odds = c("$\\underline{P(Yi=flew\\,twice)}$", "$P(Yi=did\\, not\\, fly)$", " ", " ")
key3 = getting_odds(Odds, effect, pars, MASS_PER_ML[[2]], 3) 

# -1 vs. 0
#Odds = c("$\\frac{P(Yi=flew\\, in\\, T1\\, only)}{P(Yi=did\\, not\\, fly)}$", " ", " ", " ")
Odds = c("$\\underline{P(Yi=flew\\, in\\, T1\\, only)}$", "$P(Yi=did\\, not\\, fly)$", " ", " ")
key4 = getting_odds(Odds, effect, pars, MASS_PER_ML[[2]], 1) 

# 1 vs. 0
#Odds = c("$\\frac{P(Yi=flew\\, in\\, T2\\, only)}{P(Yi=did\\, not\\, fly)}$", " ", " ", " ")
Odds = c("$\\underline{P(Yi=flew\\, in\\, T2\\, only)}$", "$P(Yi=did\\, not\\, fly)$", " ", " ")
key5 = getting_odds(Odds, effect, pars, MASS_PER_ML[[2]], 2) 

# -1 vs. 1
#Odds = c("$\\frac{P(Yi=flew\\, in\\, T1\\, only)}{P(Yi=flew\\, in\\, T2\\, only)}$", " ", " ", " ")
Odds = c("$\\underline{P(Yi=flew\\, in\\, T1\\, only)}$", "$P(Yi=flew\\, in\\, T2\\, only)$", " ", " ")
key6 = getting_odds(Odds, effect, pars, MASS_PER_ML[[3]], 2)

key = rbind(key1, key2, key3, key4, key5, key6)
rownames(key)<-NULL
colnames(key) = c("Odds","Effect", "Parameter", "Estimate*", "SE", "exp($\\beta$)*", "Wald", "p")
kable(key, caption="Table 1. Estimated Parameters, Standard Errors, and Wald Test Statistics for All Main Effects Model") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))  %>%
  column_spec(1, width="6cm", background="white") %>%
  kable_classic(html_font = "Cambria") %>%
  row_spec(c(4,8,12,16,20), extra_css = "border-bottom: 1px solid") %>%
  #column_spec(1, color = "black", background = "white") %>%
  footnote(general = " * Mass estimates were multiplied by 20 and wing-to-body estimates were divided by 100 in order to better represent how changes in the main effects would impact the odds. Instead of a 1% mass increase, which is relatively too small, or a 1 unit increase in the ratio, which is 2 magnitudes higher than approximately one standard deviation from the mean. ",
           general_title = " ", number_title = "* "
           )
Table 1. Estimated Parameters, Standard Errors, and Wald Test Statistics for All Main Effects Model
Odds Effect Parameter Estimate* SE exp(\(\beta\))* Wald p
\(\underline{P(Yi=flew\,twice)}\) Intercept \(\alpha\) 1.14 0.24 3.13 22.48 <0.001
\(P(Yi=flew\, in\, T1\, only)\) \(\delta\) Mass % \(\beta_1\) -0.4 0.01 0.67 7.71 0.005
Sex \(\beta_2\, (F=1 | M=-1)\) -0.19 0.2 0.83 0.88 0.348
Wing-to-Body \(\beta_3\,(\mu=0)\) 0.04 10.81 1.04 0.15 0.696
\(\underline{P(Yi=flew\,twice)}\) Intercept \(\alpha\) 7.18 0.2 1312.91 1340.92 <0.001
\(P(Yi=flew\, in\, T2\, only)\) \(\delta\) Mass % \(\beta_1\) 0.4 0.02 1.49 0.91 0.34
Sex \(\beta_2\, (F=1 | M=-1)\) 5 0.2 148.41 635.31 <0.001
Wing-to-Body \(\beta_3\,(\mu=0)\) 0.33 18.76 1.39 3.07 0.08
\(\underline{P(Yi=flew\,twice)}\) Intercept \(\alpha\) 0.2 0.17 1.22 1.38 0.241
\(P(Yi=did\, not\, fly)\) \(\delta\) Mass % \(\beta_1\) 0.4 0.01 1.49 4.59 0.032
Sex \(\beta_2\, (F=1 | M=-1)\) -0.76 0.17 0.47 21.07 <0.001
Wing-to-Body \(\beta_3\,(\mu=0)\) 0.28 9.72 1.32 8.36 0.004
\(\underline{P(Yi=flew\, in\, T1\, only)}\) Intercept \(\alpha\) -0.94 0.24 0.39 14.85 <0.001
\(P(Yi=did\, not\, fly)\) \(\delta\) Mass % \(\beta_1\) 0.8 0.01 2.23 18.1 <0.001
Sex \(\beta_2\, (F=1 | M=-1)\) -0.57 0.21 0.57 7.28 0.007
Wing-to-Body \(\beta_3\,(\mu=0)\) 0.24 12.06 1.27 3.88 0.049
\(\underline{P(Yi=flew\, in\, T2\, only)}\) Intercept \(\alpha\) -8.18 0.19 0 1915.51 <0.001
\(P(Yi=did\, not\, fly)\) \(\delta\) Mass % \(\beta_1\) -0.2 0.03 0.82 0.05 0.83
Sex \(\beta_2\, (F=1 | M=-1)\) -6.95 0.19 0 1376.55 <0.001
Wing-to-Body \(\beta_3\,(\mu=0)\) -0.07 18.79 0.93 0.12 0.726
\(\underline{P(Yi=flew\, in\, T1\, only)}\) Intercept \(\alpha\) 6.05 0.23 424.11 663.4 <0.001
\(P(Yi=flew\, in\, T2\, only)\) \(\delta\) Mass % \(\beta_1\) 1 0.03 2.72 3.38 0.066
Sex \(\beta_2\, (F=1 | M=-1)\) 5.19 0.22 179.47 558.52 <0.001
Wing-to-Body \(\beta_3\,(\mu=0)\) 0.29 20.22 1.34 2.02 0.155
* Mass estimates were multiplied by 20 and wing-to-body estimates were divided by 100 in order to better represent how changes in the main effects would impact the odds. Instead of a 1% mass increase, which is relatively too small, or a 1 unit increase in the ratio, which is 2 magnitudes higher than approximately one standard deviation from the mean. 
# https://haozhu233.github.io/kableExtra/awesome_table_in_html.html
Estimated Parameters, Standard Errors, and Wald Test Statistics for All Main Effects Model
Odds Effect Parameter Estimate* SE exp(\(\beta\))* Wald p
\(\underline{P(Yi=flew\,twice)}\) Intercept \(\alpha\) 1.14 0.24 3.13 22.48 <0.001
\(P(Yi=flew\, in\, T1\, only)\) \(\delta\) Mass % \(\beta_1\) -0.4 0.01 0.67 7.71 0.005
Sex \(\beta_2\, (F=1 | M=-1)\) -0.19 0.2 0.83 0.88 0.348
Wing-to-Body \(\beta_3\,(\mu=0)\) 0.04 10.81 1.04 0.15 0.696
\(\underline{P(Yi=flew\,twice)}\) Intercept \(\alpha\) 7.18 0.2 1312.91 1340.92 <0.001
\(P(Yi=flew\, in\, T2\, only)\) \(\delta\) Mass % \(\beta_1\) 0.4 0.02 1.49 0.91 0.34
Sex \(\beta_2\, (F=1 | M=-1)\) 5 0.2 148.41 635.31 <0.001
Wing-to-Body \(\beta_3\,(\mu=0)\) 0.33 18.76 1.39 3.07 0.08
\(\underline{P(Yi=flew\,twice)}\) Intercept \(\alpha\) 0.2 0.17 1.22 1.38 0.241
\(P(Yi=did\, not\, fly)\) \(\delta\) Mass % \(\beta_1\) 0.4 0.01 1.49 4.59 0.032
Sex \(\beta_2\, (F=1 | M=-1)\) -0.76 0.17 0.47 21.07 <0.001
Wing-to-Body \(\beta_3\,(\mu=0)\) 0.28 9.72 1.32 8.36 0.004
\(\underline{P(Yi=flew\, in\, T1\, only)}\) Intercept \(\alpha\) -0.94 0.24 0.39 14.85 <0.001
\(P(Yi=did\, not\, fly)\) \(\delta\) Mass % \(\beta_1\) 0.8 0.01 2.23 18.1 <0.001
Sex \(\beta_2\, (F=1 | M=-1)\) -0.57 0.21 0.57 7.28 0.007
Wing-to-Body \(\beta_3\,(\mu=0)\) 0.24 12.06 1.27 3.88 0.049
\(\underline{P(Yi=flew\, in\, T2\, only)}\) Intercept \(\alpha\) -8.18 0.19 0 1915.51 <0.001
\(P(Yi=did\, not\, fly)\) \(\delta\) Mass % \(\beta_1\) -0.2 0.03 0.82 0.05 0.83
Sex \(\beta_2\, (F=1 | M=-1)\) -6.95 0.19 0 1376.55 <0.001
Wing-to-Body \(\beta_3\,(\mu=0)\) -0.07 18.79 0.93 0.12 0.726
\(\underline{P(Yi=flew\, in\, T1\, only)}\) Intercept \(\alpha\) 6.05 0.23 424.11 663.4 <0.001
\(P(Yi=flew\, in\, T2\, only)\) \(\delta\) Mass % \(\beta_1\) 1 0.03 2.72 3.38 0.066
Sex \(\beta_2\, (F=1 | M=-1)\) 5.19 0.22 179.47 558.52 <0.001
Wing-to-Body \(\beta_3\,(\mu=0)\) 0.29 20.22 1.34 2.02 0.155

Need to get the intercept SE and values

Odds

exp(0.02*20 - 0.76 + 28.09/100) # 20% increase in mass, female, and 1/20 increase in the wing2body ratio length
## [1] 0.9239475
exp(0.02*20 + 0.76 + 28.09/100)
## [1] 4.224496
# confidence interval for wing2body ratio
exp((28.09 + 1.96*9.72)/100) 
## [1] 1.602255
exp((28.09 - 1.96*9.72)/100)
## [1] 1.094599

Predicted Probabilities

head(pp <- fitted(model))
##           0         -1            1         2
## 1 0.7470322 0.03581826 2.459149e-07 0.2171493
## 2 0.8116776 0.02845344 2.925412e-07 0.1598686
## 3 0.6983200 0.04316854 2.166785e-07 0.2585113
## 4 0.8013507 0.03101678 2.841098e-07 0.1676322
## 5 0.7663040 0.04010556 2.581516e-07 0.1935902
## 6 0.7890842 0.03734506 2.743930e-07 0.1735705

Amazing to see what a clear impact wing2body ratio has for whether bugs are more likely to fly twice or not fly at all. Those with wing2body ratios above the mean wing2body ratio appear to be more likely to fly twice and less likely to not fly. This is true across females and males.

Females only

Baseline

df <- d %>%
  filter(!is.na(egg_diff), !is.na(mass_diff), !is.na(flew_diff), sex_c == 1)

df <- df[with(df, order(mass_diff)),]

n_tfemales = nrow(df)

df$flight_case <- relevel(as.factor(df$flight_case), ref = "0")

Barplot

Egg Case

Delta Egg Response Key
Event Encoding
laid eggs in both trials 2
laid eggs in T2 only 1
laid eggs in neither trials 0
laid eggs in T1 only -1

Baseline

df <- d %>%
  filter(!is.na(egg_case), !is.na(mass_diff), !is.na(flew_diff), sex_c == 1)

df <- df[with(df, order(mass_diff)),]

n_tfemales = nrow(df)

df$flight_case <- relevel(as.factor(df$flight_case), ref = "0")

Null model

df$flight_case = droplevels(df$flight_case) # no female bug only flew in T2
null <- multinom(flight_case ~ 1, data = df) 
## # weights:  6 (2 variable)
## initial  value 102.170943 
## final  value 93.055466 
## converged

Multinom model: flew_diff ~ egg_case

model <- multinom(flight_case ~ egg_case, data = df) 
## # weights:  9 (4 variable)
## initial  value 102.170943 
## final  value 84.114906 
## converged
model_table = calculate_P(model)
## Warning in cbind(s$coefficients, c(s$edf, s$edf, s$edf), s$standard.errors[, :
## number of rows of result is not a multiple of vector length (arg 2)
## 
##  AIC:  176.2298 
##    (Intercept)  Estimate DF Std. Err.      z   wald P > |z|
## -1      -0.194    -0.656  4     0.345 -1.902  3.616   0.057
## 2        0.697    -1.189  4     0.314 -3.782 14.306   0.000
anova(null, model, test="Chisq") # adding egg_case improves fit
## Likelihood ratio tests of Multinomial Models
## 
## Response: flight_case
##      Model Resid. df Resid. Dev   Test    Df LR stat.      Pr(Chi)
## 1        1       184   186.1109                                   
## 2 egg_case       182   168.2298 1 vs 2     2 17.88112 0.0001309677

Significant models

get_significant_modelsf(7) #-1/0 egg case not sig | 0/-1 or 2/-1 egg case not sig | -1/2 egg case not sig | but 0/2 and 2/0 egg case**
## Warning in cbind(s$coefficients, c(s$edf, s$edf, s$edf), s$standard.errors[, :
## number of rows of result is not a multiple of vector length (arg 2)

## Warning in cbind(s$coefficients, c(s$edf, s$edf, s$edf), s$standard.errors[, :
## number of rows of result is not a multiple of vector length (arg 2)

## Warning in cbind(s$coefficients, c(s$edf, s$edf, s$edf), s$standard.errors[, :
## number of rows of result is not a multiple of vector length (arg 2)

Comparing models: adding host

Adding host with egg_diff and mass_diff

Host Plant Key
Host Encoding
Golden Rain Tree (GRT) 1
Balloon Vine (BV) -1 
data <- data.frame(R = df$flight_case, 
         A = df$egg_case,
         B = df$mass_diff,
         C = df$host_c)
source("src/compare_models.R")
model_comparisonsAIC("src/generic multinomial models- multinom 1RF + 3 FF.R")
##        [,1]      [,2]      [,3]      [,4]       [,5]      
## AICs   165.1826  166.6029  167.0497  167.9211   167.9513  
## models 7         13        4         11         16        
## probs  0.3572256 0.1756085 0.1404487 0.09084268 0.08948043
## 
## m7   multinom(formula = R ~ A + B + C, data = data, trace = FALSE)
## m13  multinom(formula = R ~ B * C + A, data = data, trace = FALSE)
## m4   multinom(formula = R ~ A + B, data = data, trace = FALSE)
## m11  multinom(formula = R ~ A * B + C, data = data, trace = FALSE)
## m16  multinom(formula = R ~ B * C + A * B, data = data, trace = FALSE)
anova(m4, m7, test="Chisq") # Adding C does not improve fit
## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##       Model Resid. df Resid. Dev   Test    Df LR stat.    Pr(Chi)
## 1     A + B       180   155.0497                                 
## 2 A + B + C       178   149.1826 1 vs 2     2  5.86705 0.05320915
anova(m7, m13, test="Chisq") # Adding  mass_diff*host does not improve fit
## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##       Model Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1 A + B + C       178   149.1826                                
## 2 B * C + A       176   146.6029 1 vs 2     2 2.579779 0.2753012
model <- multinom(flight_case ~ mass_diff + egg_case, data = df) 
## # weights:  12 (6 variable)
## initial  value 102.170943 
## iter  10 value 81.763896
## iter  20 value 77.915648
## iter  30 value 77.575673
## iter  40 value 77.528423
## iter  50 value 77.524883
## final  value 77.524846 
## converged
model_table = calculate_P2(model, "mass_diff", "egg_case")
## Warning in cbind(s$coefficients, c(s$edf, s$edf, s$edf), s$standard.errors[, :
## number of rows of result is not a multiple of vector length (arg 2)
## 
##  AIC:  167.0497 
##    (Intercept) mass_diff egg_case DF    SE1   SE2    z1     z2 ML wald sex wald
## -1      -0.879    57.434   -0.534  6 17.992 0.384 3.192 -1.391  10.190    1.935
## 2        0.527    18.699   -1.089  6 14.416 0.296 1.297 -3.674   1.682   13.500
##    P1 > |z| P2 > |z|
## -1    0.001    0.164
## 2     0.195    0.000

Significant eqs

run_multinom_model = function(d) {
  m <- multinom(flight_case ~ mass_diff + egg_case, trace=FALSE, data = d) 
  model_table = calculate_P2(m, "mass_diff", "egg_case", print_table=FALSE)
  return(model_table)
}
get_significant_modelsf(11) # mass diff
## Warning in cbind(s$coefficients, c(s$edf, s$edf, s$edf), s$standard.errors[, :
## number of rows of result is not a multiple of vector length (arg 2)

## Warning in cbind(s$coefficients, c(s$edf, s$edf, s$edf), s$standard.errors[, :
## number of rows of result is not a multiple of vector length (arg 2)

## Warning in cbind(s$coefficients, c(s$edf, s$edf, s$edf), s$standard.errors[, :
## number of rows of result is not a multiple of vector length (arg 2)

get_significant_modelsf(12) # egg_case
## Warning in cbind(s$coefficients, c(s$edf, s$edf, s$edf), s$standard.errors[, :
## number of rows of result is not a multiple of vector length (arg 2)

## Warning in cbind(s$coefficients, c(s$edf, s$edf, s$edf), s$standard.errors[, :
## number of rows of result is not a multiple of vector length (arg 2)

## Warning in cbind(s$coefficients, c(s$edf, s$edf, s$edf), s$standard.errors[, :
## number of rows of result is not a multiple of vector length (arg 2)

tb = model_table
# -1 / 2 | Flew in T1, not T2
I = (tb[1,1] - tb[2,1])
M = (tb[1,2] - tb[2,2])
E = (tb[1,3] - tb[2,3])
EQ1 = paste0("log(pi_-1 / pi_2) = ", round(I, 2), " + ", round(M,2), " Mass Diff", " + ", round(E, 2), " Egg Case", "     Flew in T1, rather than both" )
EQ1  # mass_dff** | sex not
## [1] "log(pi_-1 / pi_2) = -1.41 + 38.73 Mass Diff + 0.56 Egg Case     Flew in T1, rather than both"

Predicted Probabilities

head(pp <- fitted(model))
##           0          -1          2
## 1 0.9145007 0.009854788 0.07564450
## 2 0.2856071 0.021532897 0.69286000
## 3 0.7647877 0.021005983 0.21420635
## 4 0.7482798 0.025860598 0.22585957
## 5 0.8863560 0.020152630 0.09349132
## 6 0.8810560 0.022470525 0.09647351

matrix = rbind(cbind(df$mass_diff[eggT1_rows], 
            pp[eggT1_rows,1],pp[eggT1_rows,2], pp[eggT1_rows,3]),
      cbind(df$mass_diff[egg_0rows], pp[egg_0rows,1], pp[egg_0rows,2],
            pp[egg_0rows,3]), 
      cbind(df$mass_diff[egg_2rows], pp[egg_2rows,1], pp[egg_2rows,2],
            pp[egg_2rows,3]),
      cbind(df$mass_diff[eggT2_rows], pp[eggT2_rows,1], pp[eggT2_rows,2], 
            pp[eggT2_rows,3]))
colnames(matrix) = c("mass_diff", "No", "T1", "Twice")
#write.csv(matrix, "data/animate_eggs.csv") # for interactive graphs
  • Eggs laid twice: most likely to not fly unless gain more than about 0.025 g mass and then fly only in T1 (largest mass change)
  • Eggs laid in T2: most likely to not fly (2nd widest mass change)
  • Eggs laid only in T1: most likely to fly twice (mostly only lost mass)
  • Eggs not laid: most likely to fly twice always (smallest mass change and mostly only gained about 0.030 g)

wing2body

Baseline

df <- d %>%
  filter(!is.na(egg_case), !is.na(mass_diff), !is.na(flew_diff), !is.na(wing2body), sex_c == 1)

df <- df[with(df, order(mass_diff)),]

n_tfemales = nrow(df)

df$flight_case <- relevel(as.factor(df$flight_case), ref = "0")
data <- data.frame(R = df$flight_case, 
         A = df$egg_case,
         B = df$mass_diff,
         C = df$wing2body)
source("src/compare_models.R")
model_comparisonsAIC("src/generic multinomial models- multinom 1RF + 3 FF.R") # strange error...
##        [,1]      [,2]      [,3]      [,4]      [,5]       [,6]      
## AICs   165.9791  167.0497  168.793   169.6042  169.9652   169.9733  
## models 7         4         12        11        13         8         
## probs  0.4170701 0.2441983 0.1021376 0.0680834 0.05683876 0.05661042
## 
## m7   multinom(formula = R ~ A + B + C, data = data, trace = FALSE)
## m4   multinom(formula = R ~ A + B, data = data, trace = FALSE)
## m12  multinom(formula = R ~ A * C + B, data = data, trace = FALSE)
## m11  multinom(formula = R ~ A * B + C, data = data, trace = FALSE)
## m13  multinom(formula = R ~ B * C + A, data = data, trace = FALSE)
## m8   multinom(formula = R ~ A * B, data = data, trace = FALSE)
anova(m4, m7, test="Chisq") # adding wing2body does not include fit
## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##       Model Resid. df Resid. Dev   Test    Df LR stat.    Pr(Chi)
## 1     A + B       180   155.0497                                 
## 2 A + B + C       178   149.9791 1 vs 2     2 5.070548 0.07924002
anova(m7, m12, test="Chisq") # Adding A*C does not improve fit
## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##       Model Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1 A + B + C       178   149.9791                                
## 2 A * C + B       176   148.7930 1 vs 2     2 1.186134 0.5526299
anova(m7, m11, test="Chisq") 
## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##       Model Resid. df Resid. Dev   Test    Df  LR stat.   Pr(Chi)
## 1 A + B + C       178   149.9791                                 
## 2 A * B + C       176   149.6042 1 vs 2     2 0.3749583 0.8290464
anova(m7, m13, test="Chisq") 
## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##       Model Resid. df Resid. Dev   Test    Df   LR stat.   Pr(Chi)
## 1 A + B + C       178   149.9791                                  
## 2 B * C + A       176   149.9652 1 vs 2     2 0.01392852 0.9930599
anova(m4, m8, test="Chisq") # Adding A*B does not improve fit
## Likelihood ratio tests of Multinomial Models
## 
## Response: R
##   Model Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
## 1 A + B       180   155.0497                                
## 2 A * B       178   153.9733 1 vs 2     2 1.076425 0.5837908

wing2body is not in the top model for females

Encodings

Delta Flight Response Key
Event Encoding
flew in both trials 2
flew in T2 only 1
flew in neither trials 0
flew in T1 only -1
Delta Mass Key
Event Sign
gained mass from T1 to T2 pos
no mass change between trails 0
lost mass from T1 to T2 neg
Host Plant Key
Host Encoding
Golden Rain Tree (GRT) 1
Balloon Vine (BV) -1

Females Only

Delta Flight Response Key
Event Encoding
flew in both trials 2
flew in neither trial 0
flew in T1 only -1
Delta Egg Response Key
Event Encoding
laid eggs in both trials 2
laid eggs in T2 only 1
laid eggs in neither trials 0
laid eggs in T1 only -1

Summary

flight case ~ mass diff

flight case ~ mass per

flight case ~ mass per and sex

wing2body

females only

flight case ~ mass diff and egg case

Questions

  • Egg case and egg diff not factors. Not sure how to work around this problem.
  • Animate plots?
  • For Fall season some bugs were tested more than 2 times (up to four times), but which to pull from and from which experiments?